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Abstract 


Nomenclature 


A nonlinear modeling technique was used to 
characterize response surfaces for non-dimensional 
longitudinal aerodynamic force and moment 
coefficients, based on wind tunnel data from a 
commercial jet transport model. Data were collected 
using two experimental procedures - one based on 
modem design of experiments (MDOE), and one using 
a classical one factor at a time (OFAT) approach. The 
nonlinear modeling technique used multivariate 
orthogonal functions generated from the independent 
variable data as modeling functions in a least squares 
context to characterize the response surfaces. Model 
terms were selected automatically using a prediction 
error metric. Prediction error bounds computed from 
the modeling data alone were found to be a good 
measure of actual prediction error for prediction points 
within the inference space. Root-mean-square model 
fit error and prediction error were less than 4 percent of 
the mean response value in all cases. Efficacy and 
prediction performance of the response surface models 
identified from both MDOE and OFAT experiments 
were investigated. 
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Introduction 

A fundamental scientific activity is to attempt to 
find a mathematical description for the dependence of 
an output variable on independent variables that are 
varied during an experiment. For static tests, this 
mathematical description or model can be thought of 
geometrically as a hyper-surface, also called a response 
surface. Critical issues for successfully identifying an 
adequate response surface model from experimental 
data include the experiment design (or, how the 
independent variable values are set when measuring 
the output variable response), noise level on the 
measured output, identification of a mathematical 
model structure that can capture the functional 
dependence of the output variable on the independent 
variables, accurate estimates of unknown parameters in 
the identified model structure, and the ability of the 
identified model to predict output values for data that 
were not used to identify the response surface model. 

An important application area for experiment 
design and response surface modeling is wind tunnel 
testing. Since 1997, wind tunnel testing technology 
activities conducted and sponsored by NASA Langley 
Research Center have included applying formal 
experiment design principles to empirical 
aerodynamics 1 . These activities bring to bear on 
aerodynamic research problems the powerful 
machinery of formal experiment design first introduced 
by R.A. Fisher and associates early in the 20 th century 2 
and used successfully since then in a wide range of 
industrial, scientific, and engineering applications. 
Collectively, these methods are described at NASA 
Langley as Modem Design of Experiments (MDOE), 
after a phrase from the literature of formal experiment 
design that distinguishes these methods from what is 
commonly called classical experiment design. 

Classical experimentation methods have been 
popular for hundreds of years, and form the basis of 
conventional wind tunnel testing procedures in use 
today. The defining feature of classical testing 
methods is an error control strategy that requires each 
independent variable to be changed one at a time, 
while holding all other variables constant. This 
method, formally described in the literature of 
experiment design as One Factor At a Time (OFAT) 
testing, typically involves changing the levels of the 
independent variable under study as a monotonically 
increasing function of time. This is the basis of the 
common polar, for example, which is a popular wind 
tunnel testing data structure that consists of a series of 
angle of attack levels set sequentially in a 
monotonically increasing sequence, with all other 


independent variables (Mach number, angle of sideslip, 
etc.) held constant. The OFAT experiment design has 
been used extensively for wind tunnel tests, so that 
many wind tunnel facilities have computerized data 
collection systems designed to efficiently conduct 
experiments designed in this way. 

The focus of response surface modeling is on 
defining a relationship between the output variable 
(also called a response variable or a dependent 
variable) and the independent variables (also called 
factors) that are changed during the experiment. If this 
relationship is characterized well, it is possible to 
predict responses for any combination of independent 
variables in the range of those tested, not simply the 
specific points set in the test. Furthermore, such 
mathematical relationships, called response surface 
models, can be used to predict responses in other 
circumstances, e.g., in other tests, or in flight. The 
MDOE approach to experiment design adopts tactics 
intended to make errors independent and identically 
distributed. In this way, the effects of local single- 
point errors resulting from the combined effects of all 
inevitable failures to precisely implement desired 
experimental conditions, can be made to substantially 
cancel. This results in higher precision than can be 
achieved for single points. Furthermore, the MDOE 
approach reduces the number of data points required to 
define the response surface, so the excess data points 
are available to assess the quality of predictions made 
by the model. 

Typically, once the experimental data are 
collected, polynomials in the independent variables are 
used to model the functional dependence of the output 
variable on the independent variables, and the model 
parameters are estimated from the measured data using 
least squares linear regression 3 - 4 . Unfortunately, the 
question of which polynomial terms should be included 
in the model for a given set of data is often addressed 
by trial-and-error, or by just including all polynomial 
terms that could possibly be identified from the data, 
based on information limitations. A simple example of 
an information limitation is the general statement that a 
quadratic polynomial will fit 3 data points exactly, but 
a cubic polynomial, with its 4 model parameters, 
cannot be identified from just 3 data points. This 
problem of which modeling functions to include in the 
linear regression model, called model structure 
determination, gets more difficult as the range of the 
independent variables for the model increases or the 
complexity of the underlying functional dependency 
increases. 
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Various stepwise regression techniques 3 4 - 5 can be 
used to identify an adequate model structure, but these 
methods are iterative and require the involvement of an 
experienced analyst. Neural networks using radial 
basis functions with subspace partitioning, or back 
propagation with layered and interconnected nonlinear 
activation functions, have also been applied to the 
response surface modeling problem 6 . For this type of 
approach, there is a loss of physical insight and a 
danger of overfitting the data, because the model 
structures used contain many parameters, typically 
with no mechanism for limiting the size of the model 
other than the judgment of the analyst. 

In this work, a nonlinear multivariate orthogonal 
modeling technique 7 was used to model response 
surfaces for wind tunnel data. The technique generates 
nonlinear orthogonal modeling functions from the 
independent variable data, and uses those modeling 
functions with a predicted squared error metric to 
determine appropriate model structure. The orthogonal 
functions are generated in a manner that allows them to 
be decomposed without ambiguity into an expansion of 
ordinary multivariate polynomials. This allows the 
identified orthogonal function model to be converted to 
a multivariate ordinary polynomial expansion in the 
independent variables, which provides physical insight 
into the identified functional dependencies. 

The next section gives the problem statement and 
describes the necessary theory. Following this, the 
multivariate orthogonal function modeling method is 
applied to identify response surface models for 
non-dimensional longitudinal aerodynamic force and 
moment coefficients, based on wind tunnel data from a 
commercial jet transport model. 

Theoretical Development 

The multivariate orthogonal function modeling 
approach presented here was first developed in 
Ref. [7], 

Assume an V-dimensional vector of dependent 

variable values, y ~\y\,y’ 2 ,—,yn'\ , modeled in 
terms of a linear combination of n modeling functions 
Pj , j - \,2,...,n . Each Pj is an /V-dimensional 

vector which in general depends on the independent 
variables. Then, 


The aj, j = 1,2,..., » are constant model parameters to 

be determined, and e denotes the modeling error 
vector. Eq. (1) represents the usual mathematical 
model used to fit a response surface to measured data 
from an experiment. We put aside for the moment the 
important questions of determining how candidate 
modeling functions p } should be computed from the 

independent variables, as well as which candidate 
modeling functions should be included in Eq. (1), 
which implicitly determines n. Now define an Nxn 
matrix P, 

P =[P\’P2’-> Pn] ( 2 ) 

and let a = [a l ,a 2 ,...,a n ] T . Eq. (1) can be written as a 
standard least squares regression problem, 

y = P a + e (3) 

where y is a vector of measured dependent variable 
values and P is a matrix whose columns contain 
modeling functions of the measured independent 
variables, and a is a vector of unknown parameters. 
The variable e represents a vector of errors that are to 
be minimized in a least squares sense. The least 
squares cost function is 

J = (y - Pa) T (y - Pa)= e T £ (4) 

The parameter vector estimate that minimizes this cost 
function is computed from 1-4 

a = [/ >7 >] ’ P T y (5) 

The estimated parameter covariance matrix is 

Cov(a) = £ j(a - a)(a - a) r j = <r 2 T P j (6) 

where E is the expectation operator, and the equation 
error variance a 2 can be estimated from the residuals, 

v = y-Pa (7) 


y = a l p l +a 2 p 2 + ... + a nPn +z (1) 


cf 2 = 


(N-n) L 


(y~Pa) T (y-Pa) 


T 

V V 

(N-n) 


( 8 ) 
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and n is the number of elements in parameter vector a . 
Parameter standard errors are computed as the square 
root of the diagonal elements of the Cov{a) matrix 

A n 

from Eq. (6), using a from Eq. (8). 

Estimated model output is 

y = Pi (9) 


J^y T y-f J a 1 J (p] pj) 
j= i 


or, using Eq. (11), 


J = y T y- 



( 12 ) 


(13) 


For response surface modeling, the modeling 
functions (columns of P) are often chosen as 
polynomials in the measured independent variables. 
This approach corresponds to using the terms of a 
Taylor series expansion to approximate the functional 
dependence of the output response variable on the 
independent variables. 

If the modeling functions are instead multivariate 
orthogonal functions generated from the measured 
independent variable data, advantages accrue in the 
model structure determination for response surface 
modeling. After the model structure is determined 
using the multivariate orthogonal modeling functions, 
each retained modeling function can be decomposed 
into an expansion of ordinary polynomials in the 
independent variables. Combining like terms from this 
final step puts the final model in the standard form of a 
Taylor series expansion. It is this latter form of the 
model that provides the physical insight, particularly in 
the case of modeling non-dimensional aerodynamic 
force and moment coefficients. This is the reason that 
aircraft dynamics and control analyses are nearly 
always conducted with the assumption of this form for 
the dependence of the non-dimensional aerodynamic 
force and moment coefficients on independent 
variables such as angle of attack and Mach number. 


Eq. (13) shows that when the modeling functions 
are orthogonal, the reduction in the estimated cost 
resulting from including the term aj p } in the model 
depends only the dependent variable data y and the 
added orthogonal modeling function p< . This 

decouples the least squares modeling problem, and 
makes it possible to evaluate each orthogonal modeling 
function in terms of its ability to reduce the least 
squares model fit to the data, regardless of which other 
orthogonal modeling functions are present in the 
model. When the modeling functions Pj are instead 

polynomials in the independent variables (or any other 
non-orthogonal function set), the least squares problem 
is not decoupled, and iterative analysis is required to 
find the subset of modeling functions necessary for an 
adequate model structure. 

The orthogonal modeling functions to be 
included in the model are chosen to minimize predicted 
squared error PSE, defined by 8 


psE , (>-W) r O-W) ,, g 2 JL 

N N 


(14) 


or 


Multivariate Orthogonal Function Modeling 

Orthogonal modeling functions Pj have the 
following important property: 

P, Pj = 0 , i*j , ij =1,2,-, n (10) 

Using Eqs. (2) and (TO) in Eq. (5), the / h element of 
the estimated parameter vector a is given by 

=(pJ y)/(Pj Pj) 

Combining Eqs. (2), (4), and (9)-(10), 


PSE = — + 2 a 2 max — (15) 

N N 

where <y 2 max is the maximum variance of elements in 

the error vector e, assuming the correct model 
structure. The PSE in Eq. (15) depends on the mean 

squared fit error j/N , and a term proportional to the 
number of terms in the model, n . The latter term is a 
model overfit penalty that prevents overfitting the data 
with too many model terms, which is detrimental to 
model prediction accuracy 8 . The factor of 2 in the 
model overfit penalty term accounts for the fact that 
the PSE is being used when the model structure is not 
correct, i.e., during the model structure determination 
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stage. Ref. [8] contains further justifying statistical 
arguments and analysis for the form of PSE in 
Eqs. ( 1 4)-( 1 5). Note that while the mean squared fit 

error j/n must decrease with the addition of each 
orthogonal modeling function to the model (by Eq. 
(12) or (13)), the overfit penalty term cr 2 ax n/N 
increases with each added model term (n increases). 
Introducing the orthogonal modeling functions into the 
model in order of most effective to least effective in 
reducing the mean squared fit error (quantified by 

a 2 ( pj Pj j for the y th orthogonal modeling function) 

means that the PSE metric will always have a single 
global minimum value. Figure 1 depicts this 
graphically, using actual modeling results from one of 
the cases discussed later. Ref. [8] contains details on 
the statistical properties of the PSE metric, including 
justification for its use in modeling problems. 

For wind tunnel testing, repeated runs at the same 
test conditions are often available. If a 2 is the output 
variance estimated from measurements of the output 
for repeated runs at the same test conditions, then 

o 2 max can be estimated as 

a ’max = 25<jl o (16) 

If the output errors were Gaussian, Eq. (16) 
would correspond to conservatively placing the 
maximum output variance at 25 times the estimated 
value (corresponding to a 5a 0 maximum deviation). 
However, the output errors are in general not Gaussian, 
because PSE is used to evaluate candidate model 

structures. In addition, the estimate of a 2 may not be 
very good, because of relatively few repeated runs 
available and inevitable errors in duplicating test 
conditions. These are the reasons for choosing 5 a Q 
for the maximum deviation in Eq. (16). The model 
structure determined using PSE was found to be 

virtually the same for (J 2 max in the range: 

9a 2 o <a 2 max <\00a 2 o (17) 

Using orthogonal functions to model the 
dependent variable made it possible to evaluate the 
merit of including each modeling function individually 
as part of the model, using the predicted squared error 
PSE. Since the goal is to select a model structure with 
minimum PSE, and the PSE always has a single global 
minimum for orthogonal modeling functions, the 


model structure determination was a well-defined and 
straightforward process that could be (and was) 
automated. 

After the orthogonal modeling functions that 
minimized PSE were selected, each retained 
orthogonal function was expanded into an ordinary 
polynomial expression, and common terms in the 
ordinary polynomials were combined using double 
precision arithmetic to arrive finally at a multivariate 
model using only ordinary polynomials in the 
independent variables. This procedure is explained 
below in the Parameter Estimation section. 

Orthogonal modeling functions are useful in 
determining the model structure for the dependent 
variable using the PSE metric, by virtue of the 
properties of orthogonal modeling functions and the 
resultant decoupling of the associated least squares 
problem. The subsequent decomposition of the 
retained orthogonal functions is done to express the 
results in physically meaningful terms and to allow 
analytic differentiation for partial derivatives of the 
dependent variable with respect to the independent 
variables. 

The next section describes a procedure for using 
the independent variable data to generate orthogonal 
modeling functions, which have the orthogonality 
property of Eq. ( 1 0). 

Orthogonal Function Generation 

The technique for generating orthogonal 
functions of several independent variables based on the 
data will now be described. Each orthogonal function 
Pj in general depends on the independent variables. 

Let x j, j =1,2,..., m represent Nx 1 vectors of the m 
independent variables. Each element of the Xj 

corresponds to one data point. Assign k as a positive 
integer that serves as a label for a unique set of m 
non-negative integers . For example, if 

m= 2, it might be that &=1 corresponds to {0,0}, k= 2 
corresponds to {0,1}, £=3 corresponds to {1,0}, k = 4 
corresponds to {0,2}, k = 5 corresponds to {1,1}, and so 
on. Orthogonal function p k is associated with the 
set of m non-negative integers, and also with the 
ordinary polynomial function of the m 

independent variables: 

w k =Aj‘®A- 2 2 ®...®x r ” (18) 
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where, for example, .vp denotes an /V-dimensional 
vector with each element of the independent variable 
vector .V] raised to the r { power, and ® denotes an 
element-by-element product. Note that each w k is an 
ordinary polynomial in the independent variables. 
Define the order of w k as 


<p{k) = r l +r 2 +...+r m (19) 

The orthogonal functions of m independent 
variables are defined by the following generating 
relations: 


p ] = / (20) 

where / represents an Ns 1 vector of ones and k= 1 is 
associated with the set of m zeros, /•( = r 2 =... = r m = 0 . 

Each new set k evolves from a previous set k , related 


as follows: 




k <=> | 

r\,r 2 ,~ 

-> V-1 

' r /j ■> r p+] ’ --■> r m } 

k o | 

r \- r 2 ,~ 

•> v-1 

, fjj + 1, r^ + 1 , ...,r m j 


where // is an integer. In Eq. (21), the only difference 
between the integers in set k and those in set k is that 
the integer index for independent variable ju in set k is 

one more than in set k . By assumption, the £ th 
orthogonal function has already been generated. The 
orthogonal function p k is then generated by 

Pk = - x M ®Pk - X y) Pj (22) 

j 

with the summation over all j < k such that 
<p(k)-<p(j) < 2 . The y* are constants to be 

determined. The index k and its associated integer set 
keep track of the independent variable orders for the 
yfc th orthogonal function. Each new orthogonal function 
must be orthogonal only to the previously generated 
orthogonal functions of the same order, one order 
lower, and two orders lower to be orthogonal to the 
entire set of generated orthogonal functions. Proof of 
this was found by Weisfeld 9 . Orthogonality was 
verified numerically in the software implementation of 
the technique used in this work. 


The scalars y k are computed for each j by 

T 

multiplying both sides of Eq. (22) by Pj and invoking 
the orthogonality of the p } from Eq. (10), 



P]( X M® Pk) 

T 

Pj Pj 


(23) 


The quantities on the right in Eq. (23) are either a 
measured independent variable vector or a previously 

L 

generated orthogonal function. After the scalars y j 
are determined from Eq. (23) for each j, p k can be 
computed from Eq. (22). 

The process can be repeated to generate 
orthogonal functions of arbitrary order in the 
independent variables, subject only to limitations 
related to the information contained in the data. 
Speaking loosely, the multivariate orthogonal functions 
p k are orthogonalized versions of the corresponding 
ordinary polynomial functions w k . 


Parameter Estimation 

Once the orthogonal functions to be included in 
the model of Eq. (1) were generated by Eqs. (20), (22) 
and (23), then selected by minimizing the PSE from 
Eq. (14), each retained orthogonal function was 
decomposed into an expansion of ordinary polynomial 
functions in the independent variables. This step 
introduced negligible error, as described below. 

From the orthogonal function generating equation 
(22) and the discussion in the previous section, it can 
be deduced that for any index k, the orthogonal 
function p k is a linear combination of the h>, for all 
i < k . In other words, each orthogonal function p k 
can be expressed as a linear combination of the 
w h i = \,2,...,k , which are ordinary polynomials 
corresponding to the integer sets of all previously 
generated orthogonal functions plus the current (£ lh ) 
one. In equation form, 

Pk =i>k\*’\ + t>k2 w 2+- + b kk *>k ( 24 ) 

where the b kj , i = \,2,...,k , are constants to be 
determined. There is no question as to the model 
structure given for each p k in Eq. (24), because this 
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model structure was guaranteed by the use of Eq. (22) 
in generating each orthogonal function. 

The b ki , i- 1,2 in Eq. (24) needed for the 

decomposition of each orthogonal function were 
computed using a conventional least squares solution, 
as given in Eqs. (5) and (6). Since the form and 
number of terms needed for the ordinary polynomial 
expansion of each orthogonal function was known, the 
decomposition of the orthogonal functions into an 
ordinary polynomial expansion introduced negligible 

error - on the order of the numerical precision ('0 _,2 j 

- for each orthogonal function. It would also be 
possible to find the b kj parameters in Eq. (24) exactly, 

h 

by bookkeeping and properly combining the y j 
values. 

When all retained orthogonal functions were 
decomposed, the expansions like Eq. (24) were 
substituted into Eq. (1) for each retained orthogonal 
function, and common terms in the ordinary 
polynomials were combined using double precision 
arithmetic to arrive finally at a multivariate model 
using only ordinary’ polynomials in the independent 
variables. Ordinary polynomial terms that contributed 
less than 0.1 percent of the final model output root- 
mean-square magnitude were dropped. The computer 
program that implements the procedure described here 
can determine up to 8 th order models with a maximum 
of 10 independent variables. 

In summary, a model given by Eq. (1) of the 
dependent variable in terms of orthogonal functions 
Pj was determined using the minimum PSE criterion. 

Each p j in Eq. (I) was then decomposed into an 

expansion of ordinary polynomial functions, with the 
results combined to arrive at a model of the form: 

y - c, w, + Cj h>; +. . .+C; w. (25) 

J 'l 'l '2 H ‘lip 'n p v ■' 

with 

/|,/ 2 , e{l,2 ,.--,k max ] (26) 

where k max is the total number of orthogonal 
functions generated and are ordinary 

polynomials in the independent variables. The positive 
integer values of /],/ 2 , and n p depend on the 

particular orthogonal functions retained in the model 
structure determination stage and also on the 


subsequent decomposition of each retained orthogonal 
function in terms of ordinary polynomials. The 
number of ordinary polynomial terms n p may be 

different than n from the orthogonal function 
expansion. 

The orthogonal function model of Eq. (1) was 
useful in determining the model structure for the 
dependent variable using the minimum PSE criterion, 
by virtue of the properties of orthogonal functions and 
the resultant decoupling of the associated least squares 
problem. The subsequent decomposition of the 
retained orthogonal functions was done to express the 
results in physically meaningful terms and to allow 
analytical differentiation for derivatives of the 
dependent variable with respect to the independent 
variables. 

Experimental Data 

The wind tunnel data used in this work were 
recently acquired in the National Transonic Facility 
(NTF) at NASA Langley Research Center. The model 
was a commercial jet transport, although nothing about 
the current analytical method is limited to this class of 
aircraft, which was selected for illustrative purposes. 
Non-dimensional aerodynamic lift, drag, and pitching 
moment coefficient values were measured for various 
angles of attack and Mach numbers within the cruise 
flight envelope of the vehicle. The independent 
variables were angle of attack and Mach number. The 
response variables were non-dimensional aerodynamic 
coefficients for lift force, drag force, and pitching 
moment. 

The wind tunnel data was collected using two 
different experiment designs, MDOE and OFAT, over 
roughly the same range of angle of attack and Mach 
number. For the MDOE experiment, the independent 
variables were set according to a 2 nd order Box-Wilson 
design including far-comer points, with additional 
points from a custom 3 rd order design" that was based 
on an orthogonally blocked extension to a Box-Wilson 
design. For the OFAT experiment, angle of attack 
settings were increased monotonically with time for 
various constant Mach number settings. 

The MDOE data from the full wind tunnel test 
were divided into 9 subspaces, shown in Table 1. 
Within each subspace, independent variables were set 
according to the normalized values contained in 
Table 2. Normalized independent variable values can 
be found by mapping the independent variable values 
in engineering units for each subspace onto the interval 
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[-1,1], The normalization of each independent 
variable was implemented by 


(x -x ) 

\-*niax A mm } 

where x was the normalized value of the independent 
variable, and the independent variable range in 
engineering units was [ x min , x max ] . The inverse 
transformation was 

(x + 1) , . 

x = x min 2 ' x max ~ x min ) (-®) 

For example, a normalized value of 0.707 for 
angle of attack in subspace 2 would correspond to 
3.634 degrees in engineering units, since the angle of 
attack interval for subspace 2 was [1.5,4] deg. 
Similarly, the normalized Mach number setting -0.707 
corresponds to Mach number 0.718 for subspace 2, 
with Mach number interval [0.7,0.82]. All modeling 
and prediction in this work was done using normalized 
values of the independent variables. 

Results 

The multivariate orthogonal function nonlinear 
modeling technique was first applied to experimental 
data from the individual MDOE subspaces in Table I . 

Values for cr„ iax were found from the repeated runs at 
the normalized center points of the independent 
variable ranges, using the method described above. 
Model structure determination and parameter 
estimation was done automatically using the 
orthogonal function modeling technique, without 
iteration or analyst judgment of any kind. While the 
orthogonal function modeling software allows manual 
override by the analyst in the model structure 
determination stage, this option was not used for any of 
the results presented here. All data analysis and 
modeling was done on a 450 MHz computer running 
Matlab 5.3. The orthogonal modeling technique was 
implemented as a Matlab m-file 10 . 

Figure 2 shows the response surface model fit to 
measured lift coefficient data for the normalized 
independent variable values listed in Table 2 for 
subspace 2 of Table 1. The crosses shown in Figure 2 
are the measured data used for the modeling. The 
middle plot in Figure 3 shows the residuals for the lift 
coefficient model, which are the differences between 


the measured lift coefficient data and the response 
surface model values. The solid horizontal lines are 
the 95% confidence bounds for the prediction error of 
the response surface model, computed as two times the 
square root of the PSE from Eq. (14). The prediction 
error bound calculation therefore used modeling data 
to estimate the range for the residuals in prediction 
cases. 

The circles in Figure 3 are prediction errors 
computed for 5 data points not used for the modeling. 
This data was collected at the randomly selected 
normalized independent variable settings given in 
Table 3. For these points, the identified response 
surface model w'as used to predict the measured output. 
Figure 3 shows that the 95% confidence bound for the 
prediction error computed from PSE was a 
conservative measure of actual prediction error, 
indicated by the circles. The prediction errors were 
roughly the same magnitude as the residuals for the 
modeling data (or inference data), marked by x’s. This 
indicates that the functional dependence was captured 
well by the identified response surface model. 

The other plots in Figure 3 show similar 
information for the drag force and pitching moment 
coefficients. Table 4 contains detailed information 
about the modeling and prediction for this case. The 
first two rows of Table 4 show the mean values of the 
coefficients and estimates of the standard deviation of 
the output noise from repeated center point runs. This 
gives a measure of the output noise level and the how 
the noise magnitude compares to the mean measured 
output. Residual magnitudes close to the noise level is 
the best that any model can be expected to do. The 
third row shows the prediction error estimated from the 
modeling data, and the fourth row shows the root mean 
square of the modeling error using the independent 
variable settings in Table 2 for subspace 2 of Table 1. 
The fifth row of Table 4 shows the root mean square of 
the actual prediction error using the identified response 
surface model to predict the 5 points from Table 3. 
The prediction error estimate based on the modeling 
data was conservative (i.e., too large) by approximately 
a factor of 4 for this data. The actual prediction 
performance was nearly the best that it could be, 
considering the values for the estimated output noise 
cr 0 , compared to the actual prediction errors. The 
actual prediction errors were less than 1 percent of the 
mean values for all coefficients, indicating an excellent 
prediction capability for the identified response surface 
models. Finally, the last three rows of Table 4 indicate 
the number of orthogonal function terms in each 
identified response surface model, the number of 
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ordinary polynomial terms in each model after 
decomposing the orthogonal functions, and the 
maximum order of the ordinary polynomial terms in 
the model. 

Similar analysis was done for the other subspaces 
defined in Table 1. The results showed that prediction 
errors estimated from the modeling data were always 
conservative, but by varying amounts. Also, the model 
fit errors and prediction errors were larger for the 
higher Mach number cases (subspaces 4-9), because of 
larger output noise values, estimated by cr„ . 

Response surface modeling was repeated for 
subspace 2 using only the 2 nd order Box- Wilson and 
far-comer points from the original experiment design 
in Table 2. Results plotted in Figure 4 show that 
model fit errors and prediction errors were generally 
larger, but the prediction error estimates from 
modeling data were still conservative. The increased 
number of prediction points in Figure 4 resulted from 
treating the unused 3 rd order points in Table 2 as 
additional prediction cases. 

Figures 5 and 6 show results from response 
surface modeling based on the 2 nd order Box-Wilson 
design and the 3 rd order augmented design, excluding 
the far-comer points in Table 2. In these cases, some 
prediction errors were outside prediction error bound 
estimates based on the modeling data. All of these 
outliers were in fact far-comer points, which were 
included as prediction points for these cases. In this 
context, the far-comer points were extrapolations for 
the identified models, with prediction errors higher 
than the estimates based on the modeling data. 
Including the far-comer points in the modeling 
decreased the prediction errors for both 2 nd and 3 rd 
order designs individually. In addition, use of the 3 rd 
order experiment design with far-comer points 
decreased the prediction error compared to using the 
2 nd order Box-Wilson experiment design with far- 
comer points. When the far-comer points were not 
included, models identified from the 3 rd order 
experiment design had higher prediction errors than 
models identified from the 2 nd order Box-Wilson 
design. All of this behavior was consistent for data 
from the other subspaces in Table 1 . 

Data from a larger independent variable subspace 
was also analyzed. Specifically, all the test points in 
Table 2 for subspaces 1, 2, and 3 in Table 1 were 
combined into a single MDOE data set. Figure 7 
shows the modeling and prediction results in the same 
format as before, for all three longitudinal aerodynamic 
coefficients. Table 5 contains the modeling and 


prediction results. The prediction points were chosen 
randomly in the independent variable space 
corresponding to subspaces 1 , 2, and 3 in Table I . 

Response surface modeling was then done for the 
same independent variable space corresponding to 
subspaces 1, 2, and 3, but using data from an OFAT 
experiment design. For this experiment, several 

sweeps through the angle of attack range were made 
for 4 different Mach numbers. Table 6 contains the 
modeling and prediction results, and Figure 8 shows 
the plots for this case. The same number of data points 
(90) were used for the OFAT data set as for the MDOE 
data set. The response surface model identified from 
the OFAT data was applied to the same prediction data 
points used with the MDOE data. Comparing the 
information in Tables 5 and 6, the models identified 
from the MDOE and OFAT data were of comparable 
complexity, and the model fit errors were nearly the 
same, except for the drag force coefficient, which was 
fit better using the MDOE data. Prediction errors were 
lower for lift and drag force coefficient models 
identified from the MDOE data. The prediction error 
for the drag coefficient was significantly lower (by 
roughly a factor of 3) using MDOE data, compared to 
using OFAT data. 

Figures 7 and 8 show that the orthogonal function 
modeling technique did well with either MDOE or 
OFAT data, in terms of fitting the model data 
accurately, estimating prediction error from model 
data, and actual prediction, except for the degraded 
prediction with the drag coefficient model using OFAT 
data, noted earlier. All prediction errors and model 
errors in Tables 5 and 6 were less than 4 percent of the 
mean values for all aerodynamic coefficients. 

Finally, Table 7 contains response surface 
modeling results for subspaces 1, 2, and 3, using only 
2 nd order Box-Wilson points, including far comer 
points. These data represent a subset of the MDOE 
data used for Table 5 and Figure 7. Using the MDOE 
subset of 2 nd order points with far comers (Table 7), 
the prediction error for all output response surface 
models was increased compared to the response 
surface models identified from all MDOE data (Table 
5). This behavior was similar to what was seen in 
Figures 3 and 4 for MDOE subspace 2 data. Compared 
to the models identified from OFAT data, prediction 
errors using the models identified from the MDOE 2 nd 
order subset data were lower for CD and CL, but 
higher for CM. The number of data points used to 
identify the models in Table 7 was 60, compared to 90 
for Tables 5 and 6. 
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Concluding Remarks 

A nonlinear modeling technique based on the use 
of multivariate orthogonal modeling functions 

generated from the measured data was used to 
characterize response surfaces for wind tunnel data. A 
minimum prediction error criterion was used to select 
which multivariate orthogonal functions should be 
included in the model. When used with orthogonal 
modeling functions, this prediction error metric always 
has a single global minimum. The global minimum 
was used to automatically identify which orthogonal 
modeling functions should be included in the response 
surface model. The orthogonal functions were 

generated in a manner that allows them to be 
decomposed without ambiguity into an expansion of 
ordinary multivariate polynomials. This allowed the 
final form of the model to be converted to a 
multivariate ordinary polynomial expansion in the 
independent variables, which can provide physical 
insight into the identified functional dependencies. 

The multivariate orthogonal function modeling 
technique was applied to various subsets of wind 
tunnel data for a commercial jet transport model, and 
was shown to be capable of accurately identifying 
response surface models for experimental data with 
nonlinear effects. No iteration or analyst judgment was 
required. The identified response surface models fit 
the modeling data well in all cases. The model fit 
errors and prediction errors were smaller when the 
range of independent variables in the data was smaller 
(i.e., the inference space was smaller). Prediction 
capability of all identified models was tested. 
Prediction errors were estimated somewhat 
conservatively (i.e., overestimated) using 95% 
confidence bounds based on the prediction error metric 
used to identify the response surface models. Actual 
prediction errors were within the computed 95% 
confidence bounds, except for prediction points that 
were extrapolations outside the independent variable 
space used to identify the model. 

Response surface models were identified for a 
relatively large inference space using experimental data 
collected with both a Modem Design Of Experiments 
(MDOE) approach, and a One Factor At a Time 
(OFAT) approach. The multivariate orthogonal 
function modeling technique worked well with both 
data sets, achieving less than 4 percent error in both 
modeling and prediction for all cases. Using the same 
number of data points, drag coefficient model fit and 
predictions were significantly degraded using OFAT 
data, compared to using MDOE data. Models 
identified from MDOE data with fewer test points than 


OFAT data achieved lower prediction errors for lift 

and drag response predictions, but not for pitching 

moment. 
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Table 1 Inference Subspace Data Ranges 


Table 2 Normalized Independent Variable Values 


Inference 

Subspace 

Angle of Attack 
(deg) 

Mach Number 


min 

max 

min 

max 

1 

-2.5 

1.5 

0.7 

0.82 

2 

1.5 

4 

0.7 

0.82 

3 

4 

6 

0.7 

0.82 

4 

-2.5 

1 

0.82 

0.865 

5 

1 

3.5 

0.82 

0.865 

6 

3.5 

6 

0.82 

0.865 

7 

-2.5 

1 

0.865 

0.92 

8 

1 

3 

0.865 

0.92 

9 

3 

6 

0.865 

0.92 



Point 

Number 

Angle of 
Attack 

Mach 

Number 


1 

-0.707 

-0.707 


2 

0.707 

-0.707 


3 

-0.701 

0.707 


4 

0.707 

0.707 


5 

0 

0 

tn 

6 

0 

0 

Q 

c 

o 

7 

0 

0 

C/1 

1 

8 

0 

0 

i 

g 

CO 

9 

0 

-1 

<5 

V, 

10 

0 

1 

o 

-TD 

c 

CN 

1 1 

-1 

0 

12 

1 

0 


13 

0 

0 


14 

0 

0 


15 

0 

0 


16 

0 

0 


17 

-1 

-1 

E 

18 

1 

-1 

u 

£ 

19 

-1 

1 

Ph 

20 

1 

1 


21 

-0.559 

-0.559 


22 

0.559 

-0.559 

e 

23 

-0.559 

0.559 

i 

e 

24 

0.559 

0.559 

1) 

25 

0 

-0.791 

c 

26 

0 

0.791 

5 

27 

-0.791 

0 

m 

28 

0.791 

0 


29 

0 

0 


30 

0 

0 
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Table 3 Inference Subspace 2 
Prediction Data Points 



a (deg) 

M 

Point 

Number 

Eng. 

Units 

Norm. 

Units 

Eng. 

Units 

Norm. 

Units 

I 

1.5168 

0.9865 

0.8080 

0.8000 

2 

3.2098 

0.3679 

0.7120 

-0.8000 

3 

2.2758 

-0.3794 

0.7840 

0.4000 

4 

1.8696 

-0.7043 

0.7600 

0.0000 

5 

1.8917 

-0.6866 

0.7360 

-0.4000 


Table 4 MDOE Data, Inference Subspace 2 
Modeling and Prediction Results 



CD 

CL 

CM 

Mean Value 

0.0281 

0.4817 

-0.0678 

(xlO 4 ) 

0.414 

4.878 

1.649 

sIpse 

(x 1 0 4 ) 

1.947 

20.000 

7.362 

Model Error 
p/N (xlO 4 ) 

0.817 

6.066 

3.727 

Prediction Error 
(xlO 4 ) 

0.477 

5.140 

1.593 

Orthogonal 
Function Terms 

11 

9 

9 

Ordinary 

1 1 



Polynomial Terms 

9 

9 

Maximum Order 




of Ordinary 
Polynomial Terms 

4 

3 

3 


Table 5 MDOE Data, Subspaces 1,2, and 3 
Modeling and Prediction Results 



CD 

CL 

CM 

Mean Value 

0.032 

0.447 

-0.067 

(x 1 0 4 ) 

0.721 

9.063 

1.980 

s/PSE 

(xlO 4 ) 

3.274 

47.925 

19.936 

Model Error 
Jj/N (xlO 4 ) 

2.527 

41.274 

19.043 

Prediction Error 
(xlO 4 ) 

2.604 

46.550 

20.221 

Orthogonal 
Function Terms 

15 

13 

16 

Ordinary 

Polynomial Terms 

20 

20 

20 

Maximum Order 
of Ordinary 
Polynomial Terms 

6 

6 

6 


Table 6 OF AT Data, Subspaces 1, 2, and 3 
Modeling and Prediction Results 



CD 

CL 

CM 

Mean Value 

0.030 

0.441 

-0.069 

(xlO 4 ) 

0.721 

9.063 

1.980 

sIPSE 

(xlO 4 ) 

4.059 

49.362 

21.202 

Model Error 
p/N (xlO 4 ) 

3.494 

41.605 

20.242 

Prediction Error 
(xlO 4 ) 

7.945 

52.524 

19.942 

Orthogonal 
Function Terms 

15 

15 

17 

Ordinary 

Polynomial Terms 

17 

20 

19 

Maximum Order 
of Ordinary 
Polynomial Terms 

5 

6 

6 
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Table 7 MDOE Data, Subspaces 1 , 2, and 3 
2 nd Order Box-Wilson with Far Comers 
Modeling and Prediction Results 



CD 

CL 

CM 

Mean Value 

0.033 

0.447 

-0.067 

<*o 

(xlO 4 ) 

0.703 

9.836 

1.483 

JPSE 

(xlO 4 ) 

3.570 

48.539 

16.320 

Model Error 
p/N (xlO 4 ) 

2.667 

37.263 

15.094 

Prediction Error 
(xlO 4 ) 

3.116 

49.331 

22.664 

Orthogonal 
Function Terms 

14 

12 

21 

Ordinary 

Polynomial Terms 

19 

20 

22 

Maximum Order 




of Ordinary 

6 

6 

6 

Polynomial Terms 
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Figure 2 Response Surface Model Fit 
Lift Coefficient, Subspace 2 
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Figure 7 MDOE Data, Subspaces 1, 2, and 3 
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